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Abstract 

Self-similarity may stem from two origins: the process' increments infinite variance and/or process' 
memory. The ^-value of the Gutenberg-Richter law comes from the first origin. In the frame of natural 
time analysis of earthquake data, a fall of the ^-value observed before large earthquakes reflects an increase 
of the order parameter fluctuations upon approaching the critical point (mainshock). The increase of these 
fluctuations, however, is also influenced from the second origin of self-similarity, i.e., temporal correlations 
between earthquake magnitudes. This is supported by observations and simulations of an earthquake model. 
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A large variety of natural systems exhibit irregular and complex behavior which at first look 
seems to be erratic, but in fact possesses scale-invariant structure, for example see Refs. Ql^ |2fl. A 
stochastic process X(t) is called self- similar^] with index H > if it has the property 



X{Xt)£X H X{t) V A>0. 



(1) 



where the equality concerns the finite-dimensional distributions of the process X{t) on the right- 
and the left-hand side of the equation {not the values of the process). 

A point of crucial importance in analyzing data from complex systems that exhibit scale- 
invariant structure, is the following: In several systems this nontrivial structure stems from long- 
range temporal correlations; in other words, the self-similarity originates from the process' mem- 
ory only. This is the case for example of fractional Brownian motion. Alternatively, the self- 
similarity may solely come from the process' increments infinite variance. Such an example is 
Levy stable motion (the variance of Levy stable distributions is infinite since they have heavy 
tails[4], thus differing greatly from the Gaussian ones). In general, however, the self- similarity 
may result from both these origins[|5|], the presence of which can be in principle identified when 
analyzing the complex time series in terms of the new time domain termed natural timeQ6||. 

The evolution of seismicity is a typical example of complex time series. Several traditional 
studies were focused on the variation of the b- value of the Gutenberg-Richter (G-R) law|7|], which 
states that the (cumulative) number of earthquakes with magnitude greater than (or equal to) M, 
N(>M), occurring in a specified area and time is given by 



N(> M) = 10 



a—bM 



(2) 



where b is a constant, varying only slightly from region to region and the constant a gives the 
logarithm of the number of earthquakes with magnitude greater than zero[8]. These studies found 
that the Z?-value decreases before a large event, e.g., see Ref.[|9|] (cases where Z?-value increases 
prior to and then decreases sharply before a large event have been also reported [10]). Here, con- 
sidering that the b-value itself solely focuses on the one origin of self-similarity, and in particular 
the process' increments infinite variance, we show that, when employing natural time analysis, 
the b- value decrease before large earthquakes reflects an increase of the fluctuations of the order 
parameter of seismicity when approaching the critical point (mainshock, see below). The whole 
precursory variation of the order parameter fluctuations, however, is more complex since it cap- 
tures both origins. Temporal correlations between earthquake magnitudes also play an important 



role in this precursory variation, thus leading to more spectacular results compared to the ones 
obtained when restricting ourselves to traditional analysis of fr-value alone. 



For a time series comprising N events, we define[l 1] the natural time Xk f° r me occurrence of 
the k-th event (of energy Q^) by Xk = k/N. We then study the evolution of the pair (XkiQk) or 
(XkiPk), where p^ = Qk/Yln=\ Q» * s me normalized energy released during the k-th event. The 
quantity 3»(co) is defined by 4>(<y) = Y?k=iPkexp(iG)Xk)> where ft) stands for the natural angular 
frequency, and then evaluate the real function n(ft)) = |<E>(ft))| 2 in the low frequency limit. By 
considering the Taylor expansion n(ft)) = 1 — Kift) 2 + K2OO 4 + . . . , we find that the approach of a 
dynamical system to criticality (see Chapter 8 of Ref.lQ]) is identified by means of jq, i.e., 



N ( N 



which is the variance 



fly 



Ki = (X 2 ) ~ (X) 2 = £ PkXk - L PkXk J , (3) 

k=l \k=l 



120 of natural time weighted for pk- When Qk are independent and 



identically distributed positive random variables, we obtain the "uniform" (u) distribution of p^, 



as it was defined in Ref J13h (see also p. 122 of Ref.Ja]). In this case, all p^ vary around their mean 
value l/N (cf. since En=i Pn = 1) an d the quantity K"i results QJL3J] in K u = 1/12 for large N. 

In general, in a complex time series, in order to identify the two origins of self-similarity by 
means of natural time analysis, we focus on the expectation value £(k\) of the variance K\ of 
natural time when sliding a natural time window of length / through a time series of Qk > 0, 
k=l,2,...N. 

If self-similarity exclusively results from the process' memory, the $(k\) value should change 
to K u = 1/12 for the (randomly) shuffled data. This is the case of the Seismic Electric Signals 
(SES) activities [| 14|1 . which are series of low-frequency (< 1Hz) electric signals detected a few to 
several weeks (up to five months) before an earthquake when the stress in the focal region reaches 
a critical value (and hence long range correlations develop). For example, the three upper channels 
in FigfTJb) show three SES activities that preceded major earthquakes in southern, southwestern 
and western Greece, respectively, as depicted in the map of Fig[TJa). For the sake of comparison, 
the lowest channel shows an SES activity recorded in northern Greece (close to Thessaloniki). In 
all these four cases, the analysis of their original data lead to K\ w 0.07 (see also below), which 
turns to K u = 1/12 upon shuffling the data. On the other hand, if the self- similarity results from 
process' increments infinite variance only, ${K\) should be the same (but differing from K u ) for 
the original and the (randomly) shuffled data. Finally, when both origins of self- similarity are 
present, the relative strength of the contribution of the one origin compared to that of the other can 



be quantified on the basis of Eqs.(12) and (13) of Ref. Ill5ll (see also Ref.|(]]). 

In what remains, we focus on complex time series of seismicity. Earthquakes exhibit scaling 
relations chief among which is the aforementioned G-R law[7]. For reasons of convenience, we 
write hereafter G-R law of Eq.© into the form N(> M) oc \0~ bM . Considering that the seismic 
energy E released during an earthquake is related II 16ll to the magnitude through E oc 10 cM , where 
c is around 1.5, the latter form turns to the distribution, 



P(E)ocE-V 



(4) 



where y= l+b/1.5. Hence, b ~ 1 means that the exponent yis around y=1.6 to 1.7, see Table 2.1 
of Ref. fl]. 

The complex correlations in time, space and magnitude of earthquakes have been extensively 
studied tl 1 7l - l2 ill . The observed earthquake scaling laws [|22ll seem to indicate the existence of phe- 
nomena closely associated with the proximity of the system to a critical point (e.g., see Ref. IU8I1 
and references therein). In the frame of natural time analysis, it has been suggested lll2n (see also 
pp. 249-254 of Ref.ig]) that the order parameter of seismicity is the quantity K"i . The Kj value itself 



23] 



may lead to the determination of the occurrence time of the impending mainshockla [11 
when SES data are available. In particular, when the K\ value resulting from the natural time 
analysis of the seismicity subsequent to the SES recording becomes approximately equal to 0.070, 
the mainshock occurs within a time window of the order of one week. This has been empirically 

nnr 



observed in several cases [11111151 



23H (see also Chapter 7 of Ref. 



6fl) including the three major 



earthquakes of FigfTJa) that followed the SES activities depicted in FigJTJb). An example of the 
K"i dynamics after the recording of the SES activity depicted in the third channel of FigfTJb) until 
the occurrence of the magnitude 6.4 mainshock on June 8, 2008 (blue star in FigfTJa)) is given in 
Ref. [24]. In the lack of SES data, we have to solely rely on the fluctuations of the order parame- 
ter of seismicity. Along these lines, we investigated[25] the period before and after a significant 
mainshock. Time-series for various lengths of W earthquakes that occurred before or after the 
mainshock have been studied. The probability distribution function (pdf) P(ki) versus K"i was 
found to exhibit a bimodal feature when approaching a mainshock. To quantify this feature, we 
considered the variability of K\, which is just the ratio 

P^oW/niKi), (5) 

where o(k\) and jU(fCi) stand for the standard deviation and the mean value of Kj for sliding 
window lengths /=6-40. The bimodal feature reflects that, upon approaching the mainshock (with 



the number W of the earthquakes before mainshock decreasing), the variability of K[ should in- 
crease. This was subsequently confirmed because before the M9.0 devastating Tohoku earthquake 
in Japan on March 1 1, 201 1, the variability of K\ exhibited! 3] a dramatic increase. 

In addition, we investigated ||27|l the order parameter fluctuations, but when considering a nat- 
ural time window of a fixed-length W sliding through a seismic catalog (cf. in general the results 
of complexity measures when considering W =const complements those deduced when taking 
windows of various lengths W). For earthquakes in California and Greece, we foundj27| that when 
W becomes compatible with the lead time of the SES activities (i.e., of the order of a few months), 
the fluctuations exhibit a global minimum before the strongest mainshock that occurred during a 
25- and 10-year period, respectively. 

Let us now study the interrelation between the b- value and the variability of K\. In particular, 
we investigate the expected value of K"i when a natural time window length is sliding through 
randomly shuffled power law distributed energy bursts that obey Eq.©. In Fig|2l the pdf P(k\) 
versus K"i is plotted for several b values, an inspection of which reveals that: For high ^-values, e.g., 
for b=l.5 and 1.4, the P(JCi) versus K"i curve is almost unimodal maximizing at a value somewhat 
larger than 0.070, while for smaller b a second mode emerges close to K"i ps which reflects 
that the fluctuations of K"i are larger. The computed values of the K"i variability as a function 
of the b value are plotted in the inset of Figj2fb). The general feature of this curve is more or 
less similar to that observed for example before Tohoku earthquake[|26|]; quantitative agreement 
cannot be demanded, however, because temporal correlations between the earthquake magnitudes 
are also present which influence the observed results. This is corroborated by the following results 
obtained from the Olami-Feder-Christensen (OFC) earthquake model[28]. We preferred to employ 
this model here, since it has been studied in detail in hundreds of publications, but we clarify that 
there exist more recent ones, e.g., see Ref.fl29|] where the primary role of the fault system geometry 
is emerged. 

The OFC model runs as follows: we assign a continuous random variable zij € (0, 1) to each site 
of a square lattice, which represents the local "energy". Starting with a random initial configuration 
taken from a uniform distribution in the segment (0,1), the value Zij of all sites is simultaneously 
increased at a uniform loading rate until a site ij reaches the threshold value Zthres=^ (i-e., the 
loading Af is such that (zij) + Af = 1). This site then topples which means that Zij is reset to 
zero and an "energy" azij is passed to every nearest neighbor, where the coupling parameter a can 
take values from zero to 0.25 and is the only parameter of the model, apart from the edge length L 
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of the square lattice. If this causes a neighbor to exceed the threshold, the neighbor topples also, 
and the avalanche continues until all zu < 1- Then the uniform loading increase resumes. The 
number of topplings defines the size of an avalanche or "earthquake" and (when it is larger than 
unity k increases by one) is used as in natural time analysis. Here, we use the case of free 
boundary conditions [30] in which a varies locally (Xij = n . l +K , where rijj is the actual number 
of nearest neighbors of the site ij (for sites in the bulk nu — 4, for sites at the edges = 3 and 
for the four sites at the corners = 2) and K denotes [30J] the elastic constant of the upper leaf 
springs measured relatively to that of the other springs between blocks in the Burridge-Knopoff 
model Olll . The OFC model is obviously non-conservative for K > for which cc,- ; < 0.25 in the 
bulk (for more details on the OFC modelling see pp. 349-363 of Ref.lfjj] and references therein). 

We first study the predictability of the OFC model on the basis of the K"i variability. We con- 
sider the variability fa which is a function of the natural time index k,k=l,2,...,N = 2x 10 6 
estimated by analyzing in natural time for each k the preceding W =100 avalanches. The time 
increased probability (TIP) Q32|] (i.e., the time during which there exists a high probability for the 
occurrence of a large avalanche exceeding a given threshold) is turned on when fa > fa, where fa 
is a given threshold in the prediction. If the size Qk is greater than a target avalanche size threshold 
Q c , we have a successful prediction. For binary predictions, the prediction of events becomes a 
classification task with two types of errors: missing an event and giving a false alarm. We therefore 
chooseyj] the receiver operating characteristics (ROC) graph [|34ll to depict the prediction quality. 
This is a plot of the hit rate versus the false alarm rate, as a function of the total rate of alarms, 
which here is tuned by the threshold fa. Only if in between the hit rate exceeds the false alarm 
rate, the predictor is useful. Random predictions generate equal hit and alarm rate, and hence they 
lead to the diagonal in ROC plot. Thus, only when the points lie above this diagonal the predictor 
is useful. As an example, the ROC graphs for L = 512 and K = 1 or L = 256 and K = 2 are shown 
in Fig. [3] (the rational for choosing these two cases stems from the study of Ref.[|35j] in which it 
was shown that the OFC model with free boundary conditions exhibits in these cases -see their 
Fig. 4- avalanche size distribution that agrees with the G-R law). For every given threshold value 
fa and a target threshold Q c , we get a point in this plot, thus varying fa we get a curve. The various 
curves in Fig. |3]correspond to various values of Q c = 168, . . . , 1000 increasing from the bottom to 
the top. An inspection of this figure shows that the points in each curve lie above the diagonal and 
the excess is higher for larger values of Q c . In order to investigate the statistical validity of this 
result, we include in the same graph the results where: (a) the values of fa were randomly shuffled 
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and the shuffled predictors were used (green curves) and (b) the time-series of Qk was randomly 
shuffled and then was estimated (magenta curves); in both cases, we obtain curves which al- 
most coincide with the diagonal. This clearly demonstrates that the aforementioned excess of the 
results related with the original Qk series from the diagonal comes from the sequential order of 
avalanches and cannot be considered as chancy. 

We now proceed to the investigation of the temporal correlations between the magnitudes 
nik = logio(Qk)/l.5 obtained from the sizes Qk of the avalanches in the OFC model preceding 
a large avalanche. The results can be visualized in two examples in FigJH where we plot in blue 
the exponent a^pA of the detrended fluctuation analysis (DFA) ll36ll (along with the variability /3 
plotted in red) versus the number W of avalanches before a large avalanche (negative x semi-axis, 
x = —W). Note that DFA has already been employed in Ref.[|370 for monitoring temporal cor- 
relations before bifurcations. In the upper example, FigJUa), the value of a^FA well before the 
large avalanche, being somewhat larger than 0.5, exhibits small changes but strongly increases 
upon approaching the large avalanche, i.e., at W = 100 the value of cidfa becomes fa 0.75 which 
shows intensified temporal correlations. In the lower example, Figj4fb), well before the large 
avalanche we have cidfa ~ 0.6 showing long range temporal correlations, which first turn to anti- 
correlations upon approaching the large avalanche, e.g., aoFA ~ 0.43 at W = 400, and finally 
become random, i.e, a^FA ~ 0.5 at W = 100, just before the "mainshock". In both examples of 
FigjH the variability /3 rapidly increases upon approaching a large avalanche showing clear precur- 
sory changes in the temporal correlations between avalanches' magnitudes. A detailed statistical 
study of the OFC model (K = 1, L = 512), for W = 100, 200, . . . 1000, showed that among the 579 
large (Qk > 30, 000) avalanches, only in 30% of the cases a rapid increase of /3 upon approaching 
them is observed. This is more or less consistent with empirical observations since in Japan this 
precursory increase was observed in 8 out of 25 earthquakes (all above M7 during 1 January 1994 
to 1 1 March 201 1 with depths smaller than 700 km) !126ll . Concerning the a values, when studying 
W = 100, 200, . . . 1000, among the 579 large avalanches studied, in 76% of the cases the a value 
was found to become smaller than 0.5 (as seen in FigHJb)). 
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FIG. 1: (color online) (a) Major earthquakes in Greece on January 8, 2006 (red, magnitude M w =6.7), 
February 14, 2008 (green, M w =6.9 and 6.4) and June 8, 2008 (blue, M w =6.4) (b) Their preceding SES 
activities recorded at Pirgos (PIR) measuring station located in western Greece are shown (with the cor- 
responding color) in the upper three channels. Earthquakes with SES activities at PIR are located in the 
shaded region of (a). Furthermore, an SES activity recorded at a station in northern Greece on July 13, 
2012, is depicted in the lowest channel of (b). 




FIG. 2: (color online)The probability density function P(Ki) versus K\ for several values of b for temporally 
uncorrected events obeying Eq.©. The inset depicts the variability j8 as a function of b (the cross symbols 
refer to directly computed values, while the curve has been drawn as a guide to the eye). 
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FIG. 3: (color online) The ROC diagram for the OFC earthquake model discussed in the text: red (L = 256 
and K = 2) and blue (L = 512 and K = 1) lines. In addition, two ROC diagrams are depicted based on the 
results obtained for L = 512 and K = 1 : The green curves correspond to the case when the values of fit were 
randomly shuffled and the shuffled predictors were used, while the magenta curves when the time-series of 
Qk was randomly shuffled and then fa was estimated. 
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FIG. 4: (color online) The exponent aoFA (blue, left scale) and the variability j3 (red, right scale) versus the 
number of the avalanches preceding a large avalanche, Qk = 40,325 for (a) and 2^ = 31, 145 for (b), that 
corresponds to W = for the OFC model (K=l,L = 512). 
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